#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include "struct_imatge.h"
#include "im_ppm.h"
#include <GL/glut.h>

	// Window radius
	#define MD 16
	#define multi 4

	#define PX 66
	#define PY 231

imatge inputImage,outputImage,edgesImage, singularityImage, binaryImage, enhancedImage, skeletonisationImage, minutiaeImage;
double * orientations, * direccions;
double e[8];

void directionalMatrixFiltering(imatge input, imatge output1, imatge output2, double * orientacions )
{


/* Definition of Sobel filter in horizontal direction */
  int vertSobel[3][3] = {{ -1,  0,  1 },
						 { -2,  0,  2 },
						 { -1,  0,  1 }};
  int horSobel[3][3] = {{ 1,  2,  1 },
						{ 0,  0,  0 },
						{ -1, -2,-1 }};

  int * gradients=(int *)malloc(sizeof(int)*input->h*input->w);

  double windowSumH,windowSumV;
  int x, y, i, j;  /* Loop variable */

  printf("Initializing Memory...\n");
  for(x=0;x<input->w;x++)
	  for(y=0;y<input->h;y++)
	  {
		  output1->data[pR(y,x,output1->w)]=0;
	      output1->data[pG(y,x,output1->w)]=0;
	      output1->data[pB(y,x,output1->w)]=0;

		  output2->data[pR(y,x,output1->w)]=0;
	      output2->data[pG(y,x,output1->w)]=0;
	      output2->data[pB(y,x,output1->w)]=0;
	  }

  printf("Sobel...\n");
  for (y = 1; y < input->h - 1; y++) {
    for (x = 1; x < input->w - 1; x++) {
      windowSumH = 0.0;
      windowSumV = 0.0;
      for (j = -1; j <= 1; j++) {
	    for (i = -1; i <= 1; i++) {
	      windowSumV += vertSobel[j + 1][i + 1] * input->data[pR(y+j,x+i,input->w)];
	      windowSumH += horSobel[j + 1][i + 1] * input->data[pR(y+j,x+i,input->w)];
	    }
      }
      gradients[x+y*input->w]=sqrt(windowSumH*windowSumH+windowSumV*windowSumV);
      orientacions[x+y*input->w]=atan2(windowSumH,windowSumV);

      output2->data[pR(y,x,output1->w)]=gradients[x+y*input->w];
      output2->data[pG(y,x,output1->w)]=gradients[x+y*input->w];
      output2->data[pB(y,x,output1->w)]=gradients[x+y*input->w];
    }
  }

  // Kass and Witkin for Filtered Directional Matrix
  printf("Filtering Directional Matrix...\n");
  for (y = MD; y < input->h - MD; y++) {
    for (x = MD; x < input->w - MD; x++) {
      windowSumH = 0.0;
      windowSumV = 0.0;
      for (j = -MD; j < MD; j++) {
	    for (i = -MD; i < MD; i++) {
	      windowSumV += gradients[x+i+(y+j)*input->w]*cos(2*orientacions[x+i+(y+j)*input->w]);
	      windowSumH += gradients[x+i+(y+j)*input->w]*sin(2*orientacions[x+i+(y+j)*input->w]);
	    }
      }
      windowSumV=windowSumV/(pow(MD*4,2));
      windowSumH=windowSumH/(pow(MD*4,2));

      //unsigned char pixel;
      double atangent=atan2(windowSumH,windowSumV)/2;

      direccions[x+y*input->w]=atangent;

      int pixel=(int)(((atangent+3.14159/2.0)/3.14159)*255.0);

      //int pixelInt=128+atangent*400/(2.0*3.14159);
   	  //pixel=128+atangent*400/(2.0*3.14159);
   	  //Assuring pixel is an unsigned char
      if(pixel>255){printf("Fora de rang %d\n",pixel); pixel=255;}
      if(pixel<0){printf("Fora de rang %d\n",pixel); pixel=0;}

      //Writing pixel data in grayscale
      output1->data[pR(y,x,output1->w)]=pixel;
      output1->data[pG(y,x,output1->w)]=pixel;
      output1->data[pB(y,x,output1->w)]=pixel;
    }
  }
  //Free used memory
  free(gradients);
}

double distanciaMinima(double a, double b)
{
	double dif, absolut;

	dif = b-a;
	if (dif<0) absolut = dif*(-1);
	else absolut = dif;

	while(absolut > PI/2){
		if (dif > PI/2){
			b = b-PI;
		}else b = b+PI;
		dif = b-a;
		if (dif<0) absolut = dif*(-1);
		else absolut = dif;
	}
	return dif;
}

void singularityDetection(double * o, imatge output)
{

	int x,y,i;
	double pointcare=0;
	double d[8];
	for (y = 1; y < output->h - 1; y++) {
	    for (x = 1; x < output->w - 1; x++) {
	    	pointcare=0;
	    	d[2]=o[(y-1)*output->w+x-1];
	    	d[3]=o[(y-1)*output->w+x];
	    	d[4]=o[(y-1)*output->w+x+1];
	    	d[5]=o[(y)*output->w+x+1];
	    	d[6]=o[(y+1)*output->w+x+1];
	    	d[7]=o[(y+1)*output->w+x];
	    	d[0]=o[(y+1)*output->w+x-1];
	    	d[1]=o[(y)*output->w+x-1];

	    	int num=0;
	    	//Treiem els contorns de la imatge (aquells que al seu alrededor tinguin gra majoria de zeros).
	    	for(i=0;i<8;i++){
	    		if (d[i]<0.001) num++;
	    	}
	    	if (num>6){
	    		for(i=0;i<8;i++) d[i]=0;
	    	}


	    	for(i=0;i<8;i++)
	    	{
	    		pointcare=pointcare+distanciaMinima(d[i],d[(i+1)%8]);
	    	}

	    	if (pointcare==-PI){
	    		//DELTA
	    		output->data[pR(y,x,output->w)]=255;
	    		output->data[pG(y,x,output->w)]=0;
	    		output->data[pB(y,x,output->w)]=0;
	    	}
	    	if (pointcare==PI){
	    		//LOOP
	    		output->data[pR(y,x,output->w)]=0;
	    		output->data[pG(y,x,output->w)]=255;
	    		output->data[pB(y,x,output->w)]=0;
	    		}
	    	if ((pointcare==2*PI)||(pointcare==-2*PI)){
	    		//WHORL
	    		output->data[pR(y,x,output->w)]=0;
	    		output->data[pG(y,x,output->w)]=0;
	    		output->data[pB(y,x,output->w)]=255;
	    	}


	    }
	}
}

//Rotates x,y over cx and cy with o radiants output is xd,yd
void window(int x, int y, double o, int cx, int cy, int *xd, int *yd)
{
	*xd=cx;
	*yd=cy;

	*xd+= (x-cx) * cos(o) - (y-cy) * sin(o);
	*yd+= -(x-cx) * sin(o) - (y-cy) *cos(o);
}
//Suavitza l'histograma
void smooth(int *h)
{
	int o[16];
	int i;
	for(i=0;i<16;i++) o[i]=h[i];
	h[0]=(o[0]+o[1])/2;
	h[15]=(o[14]+o[15])/2;
	for (i=1;i<15;i++)
	{
		h[i]=(o[i-1]+o[i]+o[i+1])/4;
	}
}
//Calcula l'umbral dels pics de la imatge
int umbral(imatge i)
{
	int x,y,umbral_black=255,umbral_white=0;

	for(x=0;x<i->w;x++){
		for(y=0;y<i->h;y++){
			if(umbral_black>i->data[pR(y,x,i->w)])
				umbral_black=i->data[pR(y,x,i->w)];
			if(umbral_white<i->data[pR(y,x,i->w)])
				umbral_white=i->data[pR(y,x,i->w)];
		}
	}
	printf("Pulse Umbral is %d - %d limit is %d \n", umbral_black, umbral_white, umbral_white-(umbral_white-umbral_black)/5);
	return umbral_white-(umbral_white-umbral_black)/5;
}

//Detecta pics i els binaritza
void binary(int *h, int umbralPic)
{

	int i,npics;
	int resultat[16];
	int pujant;

	for(i=0;i<16;i++) resultat[i]=255;

	pujant=(h[0]<h[1]);

	npics=90;
	while(npics>2) // fins que quedin nomes 2 pics
	{
		smooth(h); // cada volta suavitza
		npics=0;
		for(i=1;i<15;i++)
		{
			if((!pujant)&&(h[i]<h[i+1])&&(h[i]<umbralPic))
			{
				if((i-2)>0) resultat[i-2]=0;
				if((i+2)<15) resultat[i+2]=0;
				resultat[i]=0;
				resultat[i-1]=0;
				resultat[i+1]=0;
				npics++;
			}
			pujant=(h[i]<h[i+1]);
		}

	}
	for(i=0;i<16;i++) h[i] = resultat[i];
}

void binarization(imatge input, double * o, imatge output)
{
	#define MARGE 12
	#define MF 16
	int x,y,i,j;
	int umbralPic=umbral(input);
	for(x=MARGE;x<input->w-MARGE;x++)
	{
		for(y=MARGE;y<input->h-MARGE;y++)
		{
			int histograma[MF], wx,wy;
			// X,Y es el punt a evaluar
			for(i=-MF/2;i<MF/2;i++)
			{
				histograma[i+MF/2]=0;
				for(j=-MF/2;j<MF/2;j++)
				{
					//obtenim el punt a evaluar del window desplaçat segons la orientacio
					window(x+i,y+j,o[x+y*input->w], x,y,&wx,&wy);
					histograma[i+MF/2]+=input->data[pR(wy,wx,input->w)];
				}
				histograma[i+MF/2]=histograma[i+MF/2]/MF;
			}
			//histograma completat
			binary(histograma, umbralPic);
			output->data[pR(y,x,input->w)]=histograma[MF/2];
			output->data[pG(y,x,input->w)]=histograma[MF/2];
			output->data[pB(y,x,input->w)]=histograma[MF/2];
		}
	}
}

void FillRegions(unsigned long value, unsigned long region, int x, int y, imatge i, unsigned long *map)
{
	map[x+y*i->w]=region;

		if (x>0)
			if(map[(x-1)+y*i->w] == value) FillRegions(value, region, x-1, y , i,map);
		if (y>0)
			if(map[x+(y-1)*i->w] == value) FillRegions(value, region, x, y-1 , i,map);
		if (x+1<i->w)
			if(map[(x+1)+y*i->w] == value) FillRegions(value, region, x+1, y , i,map);
		if (y+1<i->h)
			if(map[x+(y+1)*i->w] == value) FillRegions(value, region, x, y+1 , i,map);


}

void blob(imatge i, unsigned long *map)
{
	unsigned long region=1,value;

	int x,y;

	//copy the image to a larger bit zone
	for(x=0;x<i->w;x++){
		for(y=0;y<i->h;y++){
			map[x+y*i->w]=i->data[pR(y,x,i->w)];
		}
	}

	for(x=0;x<i->w;x++){
		for(y=0;y<i->h;y++){
			value=map[x+y*i->w];
			if (region<65535)
			{
				if ((value==0)||(value==255))
				{
					FillRegions(value, region, x,y,i,map);
					region++;if(region==255) region=256;
				}
			}
		}
	}
}

void histogram(unsigned long *map, imatge input, unsigned long *h)
{
	int x,y;
	unsigned long i;
	for(i=0;i<65535;i++) h[i]=0;
	for(x=0;x<input->w;x++){
		for(y=0;y<input->h;y++){
			h[map[x+y*input->w]]++;
		}
	}
}

void EraseArea(unsigned long a, unsigned long *map, imatge i, imatge ori)
{
	int x,y;
	for(x=0;x<i->w;x++){
		for(y=0;y<i->h;y++){
			if(map[x+y*i->w] == a)
			{
				if(ori->data[pR(y,x,i->w)] == 255)
				{
					i->data[pR(y,x,i->w)]=0;
					i->data[pG(y,x,i->w)]=0;
					i->data[pB(y,x,i->w)]=0;
				}
				else
				{
					i->data[pR(y,x,i->w)]=255;
					i->data[pG(y,x,i->w)]=255;
					i->data[pB(y,x,i->w)]=255;
				}
			}
		}
	}
}

void enhancement(imatge input, imatge ori)
{
	unsigned long *map;
	unsigned long h[65535],i;
	map=(unsigned long *)malloc(sizeof (unsigned long)*input->h*input->w);
	blob(ori, map);
	histogram(map, ori ,h);
	for(i=1;i<65535;i++)
	{
		if((h[i]>0)&&(h[i]<15))
			EraseArea(i,map,input,ori);
	}
	free(map);
}

void valors_pixels(imatge millora, int pos_x, int pos_y, unsigned char *p, unsigned char *n0, unsigned char *n1, unsigned char *n2, unsigned char *n3, unsigned char *n4, unsigned char *n5, unsigned char *n6, unsigned char *n7, unsigned char *negat_n1, unsigned char *negat_n3, unsigned char *negat_n5, unsigned char *negat_n7)
{
	*p = millora->data[pR(pos_y,pos_x,millora->w)]; //p és el píxel actual.
	*n1 = millora->data[pR(pos_y-1, pos_x+1, millora->w)];//n1 és píxel dret superior.
	*n2 = millora->data[pR(pos_y-1, pos_x, millora->w)];//n1 és píxel superior a p.
	*n3 = millora->data[pR(pos_y-1, pos_x-1, millora->w)];//n2 és píxel esquerra superior.
	*n4 = millora->data[pR(pos_y, pos_x-1, millora->w)];//n3 és píxel esquerra.
	*n5 = millora->data[pR(pos_y+1, pos_x-1, millora->w)];//n4 és píxel inferior esquerra.
	*n6 = millora->data[pR(pos_y+1, pos_x, millora->w)];//n5 és píxel inferior.
	*n7 = millora->data[pR(pos_y+1, pos_x+1, millora->w)];//n6 és píxel inferior dret.
	*n0 = millora->data[pR(pos_y, pos_x+1, millora->w)];//n7 és píxel dret.

	//1 és figura, 0 és fons blanc.
	if (*p == 255) *p = 0; else *p = 1;
	if (*n0 == 255) *n0 = 0; else *n0 = 1;
	if (*n1 == 255) *n1 = 0; else *n1 = 1;
	if (*n2 == 255) *n2 = 0; else *n2 = 1;
	if (*n3 == 255) *n3 = 0; else *n3 = 1;
	if (*n4 == 255) *n4 = 0; else *n4 = 1;
	if (*n5 == 255) *n5 = 0; else *n5 = 1;
	if (*n6 == 255) *n6 = 0; else *n6 = 1;
	if (*n7 == 255) *n7 = 0; else *n7 = 1;

	if (*n1 == 0) *negat_n1 = 1; else *negat_n1 = 0;
	if (*n3 == 0) *negat_n3 = 1; else *negat_n3 = 0;
	if (*n5 == 0) *negat_n5 = 1; else *negat_n5 = 0;
	if (*n7 == 0) *negat_n7 = 1; else *negat_n7 = 0;
}

void erasePixel(int pos_x, int pos_y, imatge i)
{
	i->data[pR(pos_y, pos_x, i->w)] = 255;
	i->data[pG(pos_y, pos_x, i->w)] = 255;
	i->data[pB(pos_y, pos_x, i->w)] = 255;
}

void skeletonisation(imatge millora)
{
	int x, y;
	unsigned char p, n0, n1, n2, n3, n4, n5, n6, n7, negat_n1, negat_n3, negat_n5, negat_n7;
	int B4, B0, B2, B6;
	int somethingErased=1;

	while (somethingErased)
	{
		somethingErased=0;
		//Mirar cada píxel de la imatge millorada.
		for(x=0;x<millora->w;x++){
			for(y=0;y<millora->h;y++){

				valors_pixels(millora, x, y, &p, &n0, &n1, &n2, &n3, &n4, &n5, &n6, &n7, &negat_n1, &negat_n3, &negat_n5, &negat_n7);
				//Si el píxel actual p es negre de la figura (si p == 1)...
				if(p == 1){
					//Si p és contorn esquerra... (si el seu píxel esquerra és fons).
					if (n4 == 0){
						B4 = n0*(n1+n2+n6+n7)*(n2+negat_n3)*(negat_n5+n6);
						//Si B4 es cert el·liminem p de la figura i el fiquem en blanc.
						if (B4)	{erasePixel(x, y, millora);somethingErased=1;}
					}
					//Si p és contorn dret...
					if(n0 == 0){
						B0 = n4*(n2+n3+n5+n6)*(n6+negat_n7)*(negat_n1+n2);
						if (B0)	{erasePixel(x, y, millora);somethingErased=1;}
					}
					//Si p és contorn superior...
					if(n2 == 0){
						B2 = n6*(n0+n4+n5+n7)*(n0+negat_n1)*(negat_n3+n4);
						if (B2) {erasePixel(x, y, millora);somethingErased=1;}
					}
					//Si p és contorn inferior...
					if(n6 == 0){
						B6 = n2*(n0+n1+n3+n4)*(n4+negat_n5)*(n0+negat_n7);
						if (B6) {erasePixel(x, y, millora);somethingErased=1;}
					}
				}
			}
		}
	}
}

void ValorsContorn(imatge esquelet, int pos_x, int pos_y, unsigned char *veins_p)
{
	//1 és figura, 0 és fons blanc.
	veins_p[0] = esquelet->data[pR(pos_y-1, pos_x+1, esquelet->w)];
	if (veins_p[0] == 255) veins_p[0] = 0; else veins_p[0] = 1;
	veins_p[1] = esquelet->data[pR(pos_y-1, pos_x, esquelet->w)];
	if (veins_p[1] == 255) veins_p[1] = 0; else veins_p[1] = 1;
	veins_p[2] = esquelet->data[pR(pos_y-1, pos_x-1, esquelet->w)];
	if (veins_p[2] == 255) veins_p[2] = 0; else veins_p[2] = 1;
	veins_p[3] = esquelet->data[pR(pos_y, pos_x-1, esquelet->w)];
	if (veins_p[3] == 255) veins_p[3] = 0; else veins_p[3] = 1;
	veins_p[4] = esquelet->data[pR(pos_y+1, pos_x-1, esquelet->w)];
	if (veins_p[4] == 255) veins_p[4] = 0; else veins_p[4] = 1;
	veins_p[5] = esquelet->data[pR(pos_y+1, pos_x, esquelet->w)];
	if (veins_p[5] == 255) veins_p[5] = 0; else veins_p[5] = 1;
	veins_p[6] = esquelet->data[pR(pos_y+1, pos_x+1, esquelet->w)];
	if (veins_p[6] == 255) veins_p[6] = 0; else veins_p[6] = 1;
	veins_p[7] = esquelet->data[pR(pos_y, pos_x+1, esquelet->w)];
	if (veins_p[7] == 255) veins_p[7] = 0; else veins_p[7] = 1;
}

unsigned char CalculCn(unsigned char *veins_p)
{
	unsigned char cn=0;
	int valor=0;

	for(int i = 1; i<9; i++){
		valor = veins_p[(i%8)] - veins_p[(i-1)];
		if (valor<0) valor = valor*(-1);
		cn+= valor;
	}
	cn= cn/2;
	return cn;
}

#define BIFURCATION 1
#define TERMINATION 2
/**
 * Retorna el tipus de minutia per al punt x,y
 */
int isMinutiae(imatge i,int x, int y)
{
	int r,g,b;
	r=i->data[pR(y,x,i->w)];
	g=i->data[pG(y,x,i->w)];
	b=i->data[pB(y,x,i->w)];

	if(!r && g && b) return BIFURCATION;
	if(r && !g && !b) return TERMINATION;

	return 0;
}

void minutiaeFiltering(imatge input)
{
	#define W 3
	int type,typeAdjacent,x,y,i,j;

	for(x=0;x<input->w;x++){
		for(y=0;y<input->h;y++){
			type=isMinutiae(input,x,y);
			//Check if is a minutiae
			if (type)
			{
				//Check a window arround the minutiae
				for(i=x-W;i<x+W;i++)
				{
					for(j=y-W;j<y+W;j++)
					{
						if((x==i)&&(y==j)){}
						else
						//Check bounds
						if((i>=0)&&(i<input->w)&&(j>=0)&&(j<input->h))
						{
							//determine if this is a minutiae also and is different or both are terminations
							typeAdjacent=isMinutiae(input,i,j);
							if((typeAdjacent) && ((typeAdjacent!=type) || ((typeAdjacent==type)&&(type==TERMINATION))))
							{
								//Erase both minutiae
								input->data[pR(y,x,input->w)]=0;
								input->data[pG(y,x,input->w)]=0;
								input->data[pB(y,x,input->w)]=0;
								input->data[pR(j,i,input->w)]=0;
								input->data[pG(j,i,input->w)]=0;
								input->data[pB(j,i,input->w)]=0;
							}
						}
					}
				}
			}
		}
	}
}

void minutiae(imatge esquelet)
{
	int x, y;
	unsigned char veins_p[8], cn;

	//Mirar cada píxel de la imatge esquelet.
	for(x=0;x<esquelet->w;x++){
		for(y=0;y<esquelet->h;y++){
			//Si el píxel actual és negre.
			if (esquelet->data[pR(y, x, esquelet->w)] == 0){
				ValorsContorn(esquelet, x, y, veins_p);
				cn = CalculCn(veins_p);
				//Si es Termination minutia pintar amb vermell.
				if (cn == 1){
					minutiaeImage->data[pR(y,x,minutiaeImage->w)] = 255;
					minutiaeImage->data[pG(y,x,minutiaeImage->w)] = 0;
					minutiaeImage->data[pB(y,x,minutiaeImage->w)] = 0;
				}
				//Si el píxel actual és Bifurcation minutia pintar amb blau cel
				if (cn >= 3){
					minutiaeImage->data[pR(y,x,minutiaeImage->w)] = 0;
					minutiaeImage->data[pG(y,x,minutiaeImage->w)] = 255;
					minutiaeImage->data[pB(y,x,minutiaeImage->w)] = 255;
				}
			}
		}
	}
}
/**
 * Pinta la direccional filtrada en una finestra gràfica
 */
void PintaDireccions()
{
    long x, y;
    double angle, radi;
	radi=multi/2.0;
    for (x=0;x<inputImage->w;x++)
    {
    	for(y=0;y<inputImage->h;y++)
    	{
    		angle=direccions[x+y*inputImage->w];
    		glColor3f(0.0, 0.0, 0.0);
    		glBegin(GL_LINES);
    			glVertex2f(multi*x+radi-sin(angle)*radi, (multi*inputImage->h)-multi*y+radi-cos(angle)*radi);
    			glVertex2f(multi*x+radi+sin(angle)*radi, (multi*inputImage->h)-multi*y+radi+cos(angle)*radi);
    		glEnd();
    	}
    }

}

void Pintar( void )
{
	/*1r netejem el buffer de pantalla*/
	glClear(GL_COLOR_BUFFER_BIT);

	// Pinta direccions de la imatge
	PintaDireccions();

	glutSwapBuffers();

}

void grafics()
{
	/*Inicialitzacio finestra treball -> GLUT*/
	glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGBA);		/*<- FIxeu-vos que ha canviat!*/

	/*Creem la finestra de treball i la col.loquem*/
	glutInitWindowSize(inputImage->w*multi,inputImage->h*multi);
	glutInitWindowPosition(0,0);
	glutCreateWindow("Biometric Identification - Prac2");

	//Display Function
	glutDisplayFunc(Pintar);

	//initializations
	/*/Indiquem amb quin color es netejar el buffer*/
	glClearColor(1.0, 1.0, 1.0, 1.0);

	/*/Definim la matriu de projecci� 2D*/
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluOrtho2D(0.0, inputImage->w*multi, 0.0, inputImage->h*multi);

	/*/Indiquem el tamany del Viewport*/
	glViewport(0,0,inputImage->w*multi,inputImage->h*multi);

	/*/I despres deixem activada la matriu del model,
	//per a introduir-hi mes dades*/
	glMatrixMode(GL_MODELVIEW);
}

char * minutiaeAsString(int type)
{
	char * text=(char *)malloc(sizeof(char)*12);
	switch (type)
	{
		case BIFURCATION: strcpy(text, "BIFURCATION");break;
		case TERMINATION: strcpy(text, "TERMINATION");break;
	default: strcpy(text,"UNDEFINED");break;
	}
	return text;
}

#define DELTA 1
#define LOOP 2
#define WHORL 3
#define NONE 0

char * singularityAsString(int type)
{
	char * text=(char *)malloc(sizeof(char)*12);
	switch (type)
	{
		case DELTA: strcpy(text, "DELTA");break;
		case LOOP: strcpy(text, "LOOP");break;
		case WHORL: strcpy(text, "WHORL");break;
		default: strcpy(text,"UNDEFINED");break;
	}
	return text;
}

int isSingularity(imatge i, int x, int y)
{

		//DELTA
		if(		(i->data[pR(y,x,i->w)]==255) &&
				(i->data[pG(y,x,i->w)]==0) &&
				(i->data[pB(y,x,i->w)]==0))
		return DELTA;


		//LOOP
		if(		(i->data[pR(y,x,i->w)]==0) &&
				(i->data[pG(y,x,i->w)]==255) &&
				(i->data[pB(y,x,i->w)]==0))
		return DELTA;

		//WHORL
		if(		(i->data[pR(y,x,i->w)]==0) &&
				(i->data[pG(y,x,i->w)]==0) &&
				(i->data[pB(y,x,i->w)]==255))
		return WHORL;

		return NONE;
}

bool modifyOrientationMinutiae(imatge i,int x, int y)
{
	//Find where is the black pixel of your contour to know from where the line comes.
	//If up from the pixel are black...
	if((i->data[pR(y-1,x,i->w)]==0)&&(i->data[pG(y-1,x,i->w)]==0)&&(i->data[pB(y-1,x,i->w)]==0)) return true;
	//Si adalt a l'esquerra del píxel hi ha línia (és negre).
	if((i->data[pR(y-1,x-1,i->w)]==0)&&(i->data[pG(y-1,x-1,i->w)]==0)&&(i->data[pB(y-1,x-1,i->w)]==0)) return true;
	//Si a l'esquerra del píxel hi ha línia (és negre).
	if((i->data[pR(y,x-1,i->w)]==0)&&(i->data[pG(y,x-1,i->w)]==0)&&(i->data[pB(y,x-1,i->w)]==0)) return true;
	//Si adalt a la dreta del píxel hi ha línia (és negre).
	if((i->data[pR(y-1,x+1,i->w)]==0)&&(i->data[pG(y-1,x+1,i->w)]==0)&&(i->data[pB(y-1,x+1,i->w)]==0)) return true;

	return false;
}

#define MARGIN 20

void writeTextOutput(char * OutputText, imatge i, imatge s, double * d)
{
	FILE * f;
	float orienta;

	f=fopen(OutputText, "w");
	for(int x=MARGIN;x<i->w-MARGIN;x++)
	{
		for(int y=MARGIN;y<i->h-MARGIN;y++)
		{
			int type=isMinutiae(i,x,y);
			if(type)
			{
				char * textType=minutiaeAsString(type);
				//If minutiae is termination, check that orientation is the correct.
				if(type==TERMINATION){
					if(modifyOrientationMinutiae(i,x,y))
					{
						orienta = d[y*i->w+x]+3.1416;
						//Orientation should be between -pi to pi.
						if (orienta > 3.1416) orienta = orienta-2*3.1416;
						if (orienta < -3.1416) orienta = orienta+2*3.1416;
						fprintf(f, "Minutiae at %d,%d - orientation %f - type %s\n",x,y,orienta,textType);
					} else fprintf(f, "Minutiae at %d,%d - orientation %f - type %s\n",x,y,d[y*i->w+x],textType);
				}else{
					fprintf(f,"Minutiae at %d,%d - orientation %f - type %s\n",x,y,d[y*i->w+x],textType);
					free(textType);
				}
			}
			type=isSingularity(s,x,y);
			if (type)
			{
				char * textType=singularityAsString(type);
				fprintf(f,"Singularity at %d,%d - orientation %f - type %s\n",x,y,d[y*i->w+x],textType);
				free(textType);
			}
		}
	}
	fclose(f);
}


int main( int argc, char** argv  )
{
	int error;
	printf("\n --Finger--  Number of arguments is: %d\n\n",argc-1);
	// Filter argument omission
	if(argc<2)
	{
		printf("Error. \nUsage: finger <input image> [-gui] (extension .ppm assumed)\n\n");
		error=1;

	}else{

		char inputFile[255], outputFiltered[255], OutputEdges[255], OutputSingular[255], OutputBinary[255], OutputEnhanced[255], OutputSkeletonisation[255], OutputMinutiae[255], OutputText[255];

		//Prepare File Names
		strcpy(inputFile,argv[1]);
		strcpy(outputFiltered,argv[1]);
		strcpy(OutputEdges,argv[1]);
		strcpy(OutputSingular,argv[1]);
		strcpy(OutputBinary,argv[1]);
		strcpy(OutputEnhanced, argv[1]);
		strcpy(OutputSkeletonisation, argv[1]);
		strcpy(OutputMinutiae, argv[1]);
		strcpy(OutputText, argv[1]);

		strcat(inputFile,".ppm");
		strcat(outputFiltered,"_filtered.ppm");
		strcat(OutputEdges,"_edges.ppm");
		strcat(OutputSingular,"_singularities.ppm");
		strcat(OutputBinary,"_binary.ppm");
		strcat(OutputEnhanced, "_enhanced.ppm");
		strcat(OutputSkeletonisation, "_skeletonisation.ppm");
		strcat(OutputMinutiae, "_minutiae.ppm");
		strcat(OutputText,"_output.txt");

		//Read input image file
		printf("Reading file: %s...",inputFile);
		inputImage=LlegirImatgePPM(inputFile,&error);
		printf("Done\n");
		if(error)
			printf("Error reading file, aborting\n");
		else{
			outputImage=CrearImatge(inputImage);
			edgesImage=CrearImatge(inputImage);
			singularityImage=CopiarImatge(inputImage);
			binaryImage=CrearImatge(inputImage);


			orientations=(double *)malloc(sizeof(double)*inputImage->h*inputImage->w);
			direccions=(double *)malloc(sizeof(double)*inputImage->h*inputImage->w);

			//filtering

			 /* Sobel filter is applied to inputImage */
			directionalMatrixFiltering(inputImage, outputImage, edgesImage, orientations );
			printf("Singularity Detection ... \n");
			singularityDetection(direccions,  singularityImage);
			printf("Binarization ... \n");
			binarization(inputImage,orientations,binaryImage);
			printf("Enhancement ... \n");
			enhancedImage=CopiarImatge(binaryImage);
			enhancement(enhancedImage, binaryImage);
			printf("Skeletonisation ...\n");
			skeletonisationImage=CopiarImatge(enhancedImage);
			skeletonisation(skeletonisationImage);
			printf("Minutiae ...\n");
			minutiaeImage=CopiarImatge(skeletonisationImage);
			minutiae(skeletonisationImage);
			printf("Filtering Minutiae ...\n");
			minutiaeFiltering(minutiaeImage);
			printf("Finished processing\n");
			//writing output
			EscriureImatgePPM(outputImage,outputFiltered,&error);
			if (error) printf("Error writing output file: %s \n",outputFiltered );
			EscriureImatgePPM(edgesImage,OutputEdges,&error);
			if (error) printf("Error writing output file: %s \n",OutputEdges );
			EscriureImatgePPM(singularityImage,OutputSingular,&error);
			if (error) printf("Error writing output file: %s \n",OutputSingular );
			EscriureImatgePPM(binaryImage,OutputBinary,&error);
			if (error) printf("Error writing output file: %s \n",OutputBinary );
			EscriureImatgePPM(enhancedImage,OutputEnhanced,&error);
			if (error) printf("Error writing output file: %s \n",OutputEnhanced );
			EscriureImatgePPM(skeletonisationImage,OutputSkeletonisation,&error);
			if (error) printf("Error writing output file: %s \n",OutputSkeletonisation );
			EscriureImatgePPM(minutiaeImage,OutputMinutiae,&error);
			if (error) printf("Error writing output file: %s \n",OutputMinutiae );
			//Write Text output
			writeTextOutput(OutputText,minutiaeImage, singularityImage,direccions);

			printf("Work completed.\nYour output is at: %s, %s, %s, %s, %s, %s, %s, %s\n\n",outputFiltered, OutputEdges, OutputSingular, OutputBinary, OutputEnhanced, OutputSkeletonisation, OutputMinutiae, OutputText);
			//init graphics
			glutInit(&argc,argv);
			grafics();
			if (strcmp(argv[2],"-gui")==0) glutMainLoop();

			//Freeing memory
			MatarImatge(inputImage);
			MatarImatge(outputImage);
			MatarImatge(edgesImage);
			MatarImatge(singularityImage);
			MatarImatge(binaryImage);
			MatarImatge(enhancedImage);
			MatarImatge(skeletonisationImage);
			MatarImatge(minutiaeImage);
		}
	}
	return error;
}
